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Abstract 

There is an ongoing debate in the hterature as to whether the effects of averaging out 
inhomogeneities ("backreaction") in Cosmology can be large enough to account for the 
acceleration of the scale factor in the FLRW models. In particular, some simple models of 
' structure formation studied in the literature seem to indicate that this is indeed possible, 

and it has also been suggested that the perturbed FLRW framework is no longer a good 



(N 



O ' approximation during structure formation, when the density contrast becomes nonlinear. In 

this work we attempt to clarify the situation to some extent, using a fully relativistic model 
^ ■ of pressureless spherical collapse. We find that whereas averaging during structure formation 

can lead to acceleration via a selective choice of averaging domains, the acceleration is not 
present when more generic domains are used for averaging. Further, we show that for most 
I of the duration of the collapse, matter velocities remain small, and the perturbed FLRW 

■ form of the metric can be explicitly recovered, in the structure formation phase. We also 

T^lj- I discuss the fact that the magnitude of the average effects of inhomogeneities depends on the 

lO ' scale of averaging, and while it may not be completely negligible on intermediate scales, it 

is expected to remain small when averaging on suitably large scales. 

o 

oo . 1 Introduction 
O 

Understanding the origin of Dark Energy is undoubtedly one of the most important problems 
^ I in cosmology today. Recently, several workers in the field (including the present authors) 

5^ ' [D El El H]) have supported the hypothesis that the acceleration of the scale factor in the 

standard Friedmann-Lemaitre-Robertson- Walker (FLRW) cosmologies, might be attributed to 
the effects of averaging over inhomogeneities in the real Universe. It is known that explicit 
averaging of the Einstein equations leads to non-trivial corrections O El [7] . These corrections 
or "backreaction", it is argued, can mimic a Dark Energy term in the Cosmological equations. 
Specifically, some simple models have been presented which demonstrate how such an effect 
may arise [1]. 

On the other hand, it has also been argued that as long as the standard perturbed FLRW 
picture for the metric of the Universe is a good approximation, any effects of averaging over the 
inhomogeneous perturbations must remain small (see, e.g. Refs. [51 [9l [lOl [11] ) . The counter 
given to this argument is that the perturbed FLRW framework is expected to break down 
around the time of structure formation [1]. That A^-body simulations do not demonstrate any 
such breakdown is attributed to the fact that these simulations always work in the Newtonian 
limit of General Relativity (GR), usually employing periodic boundary conditions, and that 
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the backreaction due to averaging under such conditions can be shown to vanish [T^]. Further, 
since the evolution of the background scale factor is fixed in the simulations, backreaction effects 
would not be visible even if present. 

This situation clearly demands some clarification. One expects that as long as gravitational 
fields are weak, the perturbed FLRW framework should work well. Is the weak field approxi- 
mation then actually breaking down during structure formation? If so, then clearly a radical 
rethink of the approach to Cosmology is in order. If not, on the other hand, then one needs 
to understand how simple models such as the one presented by Rasanen [1] can achieve an 
acceleration upon averaging during structure formation. 

In this paper we will attempt to address these issues. We employ a simple but fully rela- 
tivistic model of spherical collapse, with appropriate initial conditions imposed. We will show 
that the effect observed by Rasanen can actually be attributed to a very selective choice of 
averaging domain, and that in fact the effect is not present when a more generic averaging 
domain is chosen. However, we do find that the effect of averaging inhomogeneities may not 
be completely negligible. In particular we see ~ 10% deviations of the effective deceleration 
parameter to be defined below, from the expected FLRW value, on the scale over which we 
define our averaged quantities. This is in line with the findings of Li and Schwarz |131ll4j in the 
context of averaging of perturbative inhomogeneities. To support these results, we will show 
that at any stage of the collapse, if matter velocities are small, then the weak field conditions 
hold and one can explicitly recover the perturbed FLRW form of the metric. We will see that 
matter velocities do, in fact, remain small for most of the duration of the collapse. We will 
also present an argument explaining the origin of the ~ 10% effect from perturbative metric 
inhomogeneities on the scales at hand, and argue that in the real universe, the averaging scale 
is expected to be large enough for these effects to be much smaller. 

We have organised the paper as follows. In section 2 we will set up our model for the 
spherical collapse, and compare its results with those obtained by Rasanen. In section 3 we 
will show how the perturbed FLRW metric is recovered, and we will conclude in section 4 with 
a discussion of further tests of the backreaction argument. 

2 Spherical Collapse : Setting up the model 

Before we describe the model we use, it will help to recall the model used by Rasanen [T], which 
is what we wish to compare with. Rasanen's model can be summarized as follows : one considers 
two disjoint regions, one overdense and the other completely empty, each evolving according to 
the FLRW evolution equations. (The embedding of these regions in an FLRW background, and 
the behaviour of the region between these two regions, is not considered.) The scale factor in the 
overdense region therefore behaves as ai oc (1 — cosu) with t oc {u — sinn), and the scale factor 
in the empty region behaves as 02 oc It is then straightforward to show that if one defines a 
volume averaged scale factor by = a'^ + , then the effective deceleration parameter given by 
q = —{da) /a? becomes negative (indicating acceleration) around the time that the overdense 
region turns around and starts collapsing. In what follows, we will attempt to reconstruct such 
a situation more rigorously, with appropriate matching and initial conditions taken care of. 

2.1 The LTB solution 

We will employ a model containing spherically symmetric pressureless "dust" , with the energy- 
momentum tensor given by 



Tab = diag{p{t,r), 0,0,0) . 
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The line element of interest is hence the Lemaitre-Tolman-Bondi (LTB) metric [17] given in 
synchronous and comoving coordinates, by 

E>'2 J 2 

ds^ = -dt^ + -^^yf-^ + R'dn^ , (2) 
1 — k[r)r^ 

where R{t, r) is the area radius of a shell of matter labelled by the comoving radius r, and k{r) 
is a function of integration to be determined by the initial conditions. (A prime and a dot refer 
to derivatives with respect to r and t respectively. We set the speed of light c = 1, except when 
displaying numerical results where we will explicitly account for factors of c.) The equations of 
motion are 

where M(r) is another function of integration which can be determined from the initial density 
profile once a scaling freedom in the coordinate r is fixed. The spherical collapse model has 
been extensively studied in the literature in the context of structure formation [18]. However, 
to the best of our knowledge, the issue of whether or not the collapse situation can be recast 
as perturbed FLRW, has not been discussed. 

The solutions to Eqn. ([3|) depend on the sign of the function k{r), and can be given in 
parametric form as 



f9GM{r) 
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Rir,t) = — ^t-ts{r)f'\ forA:(r)=0, (5a) 

GMir) , ^ , , ^ GM(r) , ^ , . w n 

it = — , , , „ (cosh 77 — 1) : t — ts[r) = 5-777 (sinhjy — ri) , < 77 < 00 , for klr) < , 

(5b) 

GMir) GMir) 
^=lTTT(l-cosr?) ; t - ts{r) = ^ j (r? - sin??) , < 7? < 27r, for A;(r) > . 

(5c) 

The function ts{r) appearing above is completely determined once k{r) and M{r) are known 
and a scaling choice for r is made. Since we will always start the evolution at t = > ts{r), 
this function will be of little physical relevance, except to ensure consistency at t = ti. 



2.2 Initial conditions 

While choosing the initial density, velocity and coordinate scaling profiles, we make the impor- 
tant assumption that at initial time, a well-defined global background FLRW solution can be 
identified, with scale factor a{t), Hubble parameter H(t) and density p{t). This is reasonable 
since the CMB data (combined with the Copernican principle) assure us that inhomogeneities 
at the last scattering epoch were at the level of 10 parts per million. This assumption plays 
a crucial role in deciding which regions are overdense and will eventually collapse, and which 
regions will keep expanding (see also the first paper in Ref. [1] for a discussion on this point). 

• Initial density profile p(ti,r) : 
The initial density is chosen to be 

{(1 + 5*), r<n 
(1 - 6^), <r <r^ (6) 
1, r>r^, 
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where pi = p{ti). Initially, the region r < r^< is assumed to contain a tiny overdensity and 
the region < r < r^, an under density. In other words, 

0<5*,5^<1. (7) 

The discontinuities in the initial density profile can be smoothed out by replacing the 
step functions appropriately. We will not do this here, since the step functions make 
calculations very simple. This is not expected to affect the qualitative features of our 
final results. 

• Initial conditions on scaling and velocities : 

We match the initial velocity and coordinate scaling to the global background solution, 
by requiring 

R{ti,r) = air, (8) 
R{ti,r) = aiHir, (9) 

with and Hi denoting the initial values of the scale factor and Hubble parameter 
respectively of the global background. This amounts to setting the initial velocities to 
match the Hubble flow, ignoring initial peculiar velocities. This is only a convenient 
choice and the introduction of initial peculiar velocities is not expected to modify our 
final results qualitatively. 

For the FLRW background we consider an Einstein-deSitter (EdS) solution with scale factor 
and Hubble parameter given by 

a{t) = {t/tof/^ ; to = 2/(3i?o), (10) 
H{t) = d/a = 2/{3t), (11) 

with to denoting the present epoch. Oj fixes the initial time as 

U = 2/{3Ho)af^ . (12) 

We will always use Oj = 10~^, so that the initial conditions are being set around the CMB last 
scattering epoch; in general Oj must be treated as one of the parameters in the problem. The 
initial EdS background density is given in terms of Hq and as 

- ' Hlaf. (13) 
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2.3 Mass function M(r) and curvature function k{r) 

We now have enough information to fix M{r) and k{r). Using Eqn. (jlj) at initial time together 
with the scaling in Eqn. ([8]) gives us 



GM{r) = ]^Hlr^ { 



'1 + 5*, < r < n 

1 + {{rc/rf - l) , r, < r < r„ (14) 
^ 1 + {5Jr^) [rl -rl) , r > , 



where we have defined a "critical" radius Vc by the equation 
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The significance of will become apparent shortly. Using the initial conditions Eqns. ([5]) and 
([9]) in the evolution equation ^ at initial time, gives 

k{ry='^^M^-afHy, (16) 
air 



with Hf = H^a- ^, and hence 



k{r) = ^ { 

a. 



' S^, r < 



5^(irjrf -1) , r^<r<r^ (17) 



^ {6^/r^) (r|? -rf,) , r > r^ . 

The significance of rc is now clarified. Since 6^:,6v > 0, we have rc > r^ by definition (Eqn. 
(|15p). The following possibilities arise : 

• If Tc > r^, then k{r) > for all r, and every shell will ultimately collapse, including the 
"void" region r^ < r < r^. 

• If Tc < r„, then k{r) > for r < rc and changes sign at r = r^. Hence, the region 
r^ < r < rc will collapse even though it is underdense, while the region r > rc will expand 
forever. 

• If Tc = Tj;, then the "void" exactly compensates for the overdensity, and the universe 
is exactly EdS for r > r„. {GM{r) = (l/2)ij2^3 ^ p.] Also the "void" will 
eventually collapse. 

Clearly the most interesting case for us is the one with rc < r^j, and we will hence make this 
choice for our model. We realize that the model as it stands is not a very realistic depiction 
of the (nearly spherical) voids we see in our Universe [15], since these voids are seen to be 
surrounded by "walls" of matter. However, our goal is to describe two regions, one of which 
collapses while the other expands ever more rapidly, and our model is capable of doing so while 
retaining its fully relativistic character. 

Although we have set up the model for all values of the radial coordinate r, hereon we will 
concentrate on the region < r < r^. One reason is that most of the interesting dynamics 
takes place in this region. Another is that the region r > r„ develops shell-crossing singularities 
due to the sharp rise in density across r — r^. A more realistic model would be able to 
incorporate the pressures that are expected to build up when a shell-crossing occurs ([16j. see 
also the discussion in section 4 below), but the LTB model is limited in this respect due to its 
pressureless character. We will therefore ignore the region r > r„. 

2.4 The solution in the region < r < r„ 

The region of interest can be split into three parts : Region 1 = {0 < r < r^,}. Region 2 
= {r* < r < rc} and Region 3 = {rc < r < r„}. The solution in the three regions is as follows : 

• Region 1 (0 < r < r*) : 

R=^(^j-^r{l + 6^){l-cosu), (18a) 

tt - sin tt = ( — ) (i - ij) + (uj - sinuj) , (18b) 



1 — cosuj = - — , (18c) 



1 + (5* V Oj 
l + <5. 



R^Yi = ^ . (18d) 



r 
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For Regions 2 and 3, it is convenient to define a function e{r) as 

, 3 



1 



k{r) , < r < . 



Region 2 (r* < r < rc) : 



= i (^^^ r(l + e)(l -cosa) 



Q — sm a 



1 — cos aAr) 



1 + e \ aj 

2£ 



3/2 



(t - tj) + {ai - sinoi) , 



l + e 



re 



,3/2 



1 



e{l + e) I (1 — cosa)2 
r sma 

+ 



Hi{t — ti) sin a 



(1 + e)2 \sin ttj 



3 + e 
l + e 



(19) 

(20a) 
(20b) 
(20c) 



(20d) 



Region 3 (rc < r < r^) : 

1 / Oi 



2 V e 



r(l + e)(cosh r/ — 1) , 



sinh 7] — 7] 



2Ho [\e 



cosh?7j(r) — 1 



R^R' =^ I 



2\e\ 



3/2 



— (t - tj) + (sinh r]i - iji) , 



l + e 



re 



|3/2 



e(l+e) 



(coshr? — 1)^ 
4|e|^/^ / sinhr/ 



Hi{t — ti) sinh?7 



3 + e 
l + e 



(1 + e)^ ysinhr/j 



(21a) 

(21b) 
(21c) 



(21d) 



Tlie crossover from Region 1 to Region 2 is discontinuous in R' (but not in R) due to our 
discontinuous choice of initial density. Smoothing out the density wih also smooth out R'. The 
crossover from Region 2 to Region 3 can be shown to be smooth, by considering the limits 
r ^ r~ and r ^ or equivalently e ^ 0~ and e — > 0^. Note that the results in Eqns. (|18|) . 
(|20]) and ([2T]) are exact, and do not involve any perturbative expansions in 5* or even though 
these parameters are small. 

2.5 Behaviour of the model 

Each shell in the inner, homogeneous and overdense Region 1 behaves as a closed FLRW 
universe, expanding out to a maximum radius Rmaxi''') given by 

Rmaxir) = ^r{l + 6,) . (22) 
All the inner shells reach their maximum radius and turn around at the same time ttum given 
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Parameter name 


Parameter value 




U.UUi 


-no 


i/io.oy(jyr — (ZKm/s/ivipcj 




z/i^ojjoj — y.uD^jyr 


C 


oUD.o ivipc^jyr 


0* 


1 91^^ /"iTT //I ^2/3 9 91 V 1 n^3 


A 


yj .\J\jO 


' * 


n nn4r/ffr, — Ifi 7Mnr 


'^turni ^{j 


n 79 


Tc 


-\- Oif 1 O^j — io.o iVipC 


Tv 


1.25rc = 23.5 Mpc 




6.8 Mpc 


R{tQ,r^) 


33.3 Mpc 



Table 1: Values of various parameters used in generating plots. 



by 

Wn = t. + ^(|) (^-(n,-sm^,))-io(^) , (23) 

where we have used the smallness of and (5^, to make the last approximation. By appropriately 
choosing a value of 5^ , we can arrange for the turnaround of Region 1 to occur either before or 
after the present epoch. 

In Table [T] we have listed the parameter values which we will use frequently in displaying 
plots. Along with the parameter set {oj, ffo; (5*, t*, r^,}, we have also listed the values of the 
derived quantities {rc,ti,tQ,tturn} and speed of light c in units of MpcGyr"^. We have also 
shown the values of the present day physical area radius R{tQ,r) at r = r* and r = r^. The 
density contrasts are to be understood to reflect the inhomogeneities in the dark matter density 
close to last scattering, and not the inhomogeneities of the baryons which were much smaller 

[ia[2o]. 

In Fig. [T]we have shown the evolution of the density contrast 5{t,r) defined in the usual 
way by 

l + *«,r) = ^. (24, 

for the parameter choices of Tabled! for which one has ttum/io — 0.72, so that the collapse is 
well under way in Region 1 at the present epoch. The two panels show the contrast for two 
representative values of r, one in Region 1 and the other in Region 3. 

For clarity, we define the volume of each of our three comoving regions separately, as 

Vi= Att / =dr ; V2 = 4:1t =dr ; V3 = 4tt =dr . 

Jo yl — A;(r)r2 Jr, y 1 — A;(r)r2 Jr^ y'l — k{r)r'^ 

(25) 

The total volume of the region is used to define a volume averaged scale factor as 

<^it)^[-^y ' ^(*)=^i(^)+^2(t)+i^3(t), (26) 
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and hence an effective deceleration parameter q given by 

On the other hand, we note that Rasanen's model can be mimicked more closely by ours, if 
we simply remove the Region 2, by hand. By doing so we are left with two disjoint regions, 
each spherically symmetric, one of which is collapsing and the other expanding ever rapidly and 
becoming ever emptier. There is no physical reason to throw away Region 2 in this manner, but 
for the sake of comparison we will define a "modified" scale factor flmod and it's corresponding 
deceleration parameter Qmod by 

1 /3 

/,N _ f yi{t)+y3{t) \ _ amodamod 

Note that the normalisation of the scale factor is irrelevant in defining the deceleration pa- 
rameter. In Fig. [2] we plot q{t) and qmod{t), for several sets of initial conditions which are 
close to our "base set" listed in Tabled] (except for Fig. [2d] which has a large value for 6^j^. 
The various integrals involved in computing V{t), etc. were performed using the NIntegrate 
routine of Mathematica. To generate the plots we used the numerical derivative routine ND of 
Mathematica. The various initial conditions correspond to turnaround times that are slightly 
greater than, or slightly less than, or significantly less than the present epoch. The idea here is 
to demonstrate that the results are valid regardless of whether the collapse has just begun or is 
well under way at the present epoch. We see that while the modified scale factor does accelerate 
as in Rasanen's model, the scale factor a{t) which is the more natural choice in our situation, 
does not show this effect. The reason for this can be understood as follows. The Region 2 is 
of a rather peculiar nature - it is underdense initially and becomes emptier with time, however 
its evolution is closely linked to that of the overdense Region 1. Namely, the whole of Region 2 
(except its boundary at r = Tc), is dragged along with Region 1 and eventually turns around, 
instead of expanding away to infinity like its counterpart Region 3. Now, if one ignores Region 
2, then Rasanen's arguments about the remaining two regions stand - one region is contracting 
and the other is expanding faster than the global mean, and this stand-off leads to an acceler- 
ation of the effective scale factor a^od^ as we see in the plots of Fig. [2] But if we account for 
Region 2 as well, then we bring in a counter-balancing influence of a large underdense volume 



^Both curves in Fig. I2dl begin at g ~ 0.5 at t = ti. To enhance the contrast between the curves, we have 
plotted them for times t > O.lSttu™. The remaining plots (Figs. I2al42c|l are plotted starting from t — ti. 
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8^, = 0.001 
r. = 0.001 ( c/Hq ) 
r„ = 1.01 r. 
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1.5 

(t/t,urn) 



(b) ttur« = l.OOlSto 



8. = 1 .25 a; (371/4)'' 
8„ = 0.05 

0.005 ( c/Hq ) 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 

(t/tturn) 

(d) tturn = 0.717to 



Figure 2: The deceleration parameters for a range of parameter values. The dashed lines correspond to 
Qmod and the solid lines to q. The x-axis shows t/ttum, where ttum is the time at which Region 1 turns 
around, and is different for each plot. The values for a^, Hq and c are the same as those listed in Table 
[TJ [Both curves in Fig. [2dl begin at q ~ 0.5 at t = U.] 



which is expanding slower than average, and this reduces the accelerating influence of Region 
3 to the point of making the effect completely disappear. Note that at late times, the volume 
of Region 1 contributes negligibly to the total volume, and the volumes of regions 2 and 3 are 
comparable. 

The figures (|2ap - (pc|) show that while the full deceleration parameter q does not change sign, 
it does deviate from the EdS value of (1/2), by an amount of the order of ~ 10%. This indicates 
that while backreaction may not be large enough to cause the scale factor to accelerate, it may 
still lead to effects which are not completely negligible. 

We wish to highlight two points. First, it is very important to note the role played by 
the initial conditions in this entire excercise. The function k{r) is defined in a continuous 
fashion once the initial density, velocity and coordinate scaling are given, and k{r) then decides 
which shells will eventually collapse and which will not. The continuity of k{r) assures us that 
in models such as ours, with an overdensity surrounded by an underdensity, the underdense 
region will always contain a subregion in which k{r) > 0. We see therefore that the existence 
of Region 2, is a generic feature not restricted to our specific choice of discontinuous initial 
density or vanishing initial peculiar velocities. Further, as we see in Fig. I2dl it is possible to 
make q deviate even more significantly from the EdS value than the ~ 10% effect of the first 
three figures, but this requires an unnaturally high value of 6^ > 0.01 (the figure has 6^ = 0.05), 
which contradicts CMB data. 

Secondly, one may argue about the "naturalness" of choosing one set of regions over another 
set, in order to compute volumes. But this itself places the physicality of the acceleration effect 
into question - if one has to judiciously choose a specific set of averaging domains in order to 
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obtain acceleration on average, then the effect would appear to be an artifact of this choice 
rather than something which observers would see. 

Having seen that the average behaviour of the full region < r < r^, is close to EdS, we can 
ask whether one can explicitly show that the metric for this system is actually close to FLRW. 
We answer this in the affirmative in the next section. 



3 Transforming to Perturbed FLRW form 

We ask whether the metric ^ can be brought to the perturbed FLRW form with scalar pertur- 
bations, at any arbitrary stage of the collapse. Namely, we want a coordinate transformation 
{t, r) — > (t, f) such that the metric in the new coordinates is 

ds^ = -(1 + 2A)dP + a^ii){l + 2iP) {dr"^ + f^dQ^) , (29) 

with at least the conditions 

I i |< 1 ; I l< 1, (30) 

being satisfied. We will ignore conditions on the derivatives of ip and A for now (see the end 
of Section 3.2). The scale factor is the EdS solution, with t as the argument. The coordinate 
f is comoving with the (fictitious) background Hubble flow, but not with the matter itself. 
On physical grounds we expect that this transformation should be possible as long as the 
gravitational field is weak and matter velocity is small. We will see below that this is exactly 
what happens. In the new coordinates, all matter shells labelled by r expand with the Hubble 
flow, with a superimposed peculiar velocity. 

Since we want f to be comoving with the background, the natural choice for this coordinate 
would be f R/a, at least at early times. Also, we need to account for the local spatial 
curvature induced by the initial conditions. As an ansatz for the coordinate transformation 
therefore, we consider the equations 

f = ^^{l+at,r)), (31a) 
a{t) 

i=t + f{t,r), (31b) 
where (,{t,r) and ^^{t,r) are expected to satisfy 

I e I < 1 ; \fH\^l. (32) 

This form of the transformation keeps us close to the standard gauge transformation of cos- 
mological perturbation theory, while still accounting for the deviations in the evolution from 
the background FLRW, caused by structure formation. We will show that a self-consistent 
transformation exists, which preserves the conditions (|30|) and (j32|) for most of the evolution. 
We will use the metric transformation rule given by 

=9y(^), (33) 

and expand to leading order in the small functions ^, (,^H, A, ip and also k{r)r'^ which, as we 
see from Eqn. (jl7p . remains small in the entire region of interest. The relations in Eqn. (j33p 
must be analysed for the cases {ij) = {(tt), (tr), (rr), {69)}, in each of the three regions. (The 
remaining cases can be shown to lead to trivial or non-independent relations.) The analysis is 
similar to the standard gauge transformation analysis in relativistic perturbation theory |19j . 
Since the calculations involved are straightforward but tedious, we will only present an outline 
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of the calculation and highlight certain issues. At the end we will present equations for all three 
regions and numerically show that the transformation is well-behaved in the regime of interest. 
The cases (ij) = (99) and (rr) are easily analysed and lead to 

4, = -fH-C, (34) 

and 

e = \Kr)r' ) . (35) 

The cases (ij) = (tr) and {tt) both require | dtf |<C 1 for consistency (since the RHS of Eqn. 
(I33p in these cases has no zero order term to balance a large dtf). Note that since t is the 
proper time of each matter shell, the quantity dtf is simply the velocity of matter in the (t, f) 
frame (which is comoving with the Hubble flow). In other words, 

df 

is the radial comoving peculiar velocity of the matter shells in the (t, f) frame. We will soon 
see that whereas the quantities ^ and behave roughly as ~ {Nqv)^, the peculiar velocity av 
behaves roughly as ~ {Hqv). We will therefore treat (aw)^ as a small quantity of the same order 
as ^, etc. The case (ij) = (tr) then leads to 

^0' = avBl , (37) 

and the case {ij) = [tt) gives 

A = -i' + \{avf. (38) 



The equations (|M|) . (j35|) . (j37j) . and (|38|) are valid in the entire range < r < r„, provided the 
peculiar velocity remains small in magnitude. The comoving peculiar velocity is given by 

* = (?) • P^) 

where we have assumed for consistency that \dt{R/a)\ <^ 1 and have dropped the term {R/a)^ 
since it is expected to be of higher order than dt{R/a). (This can be seen from simple di- 
mensional considerations - we have dt{R/a) ~ HR/a, and since, from Eqns. (I35p and (llTh . 
^ ~ [HR)^, we also have {R/a)^ ~ {H R)"^ / a.) We will see that these conditions do indeed hold 
for most of the evolution, throughout the region of interest. 



3.1 The transformation in Region 1 

Since Region 1 corresponds to a homogeneous solution, the integrals in Eqns. ([35]) and (j37|) 
can be analytically performed. Eqn. (I35p directly gives 

e = \k{ry = \ (I) {H,rf , (40) 

after setting an arbitrary function of time to zero. Hence v has the structure v = ryi{t), since 
R has the structure R = ry2{t). Eqn. ([37|) then leads to 

e = haR, (41) 

after setting another arbitrary function of time to zero. [Note that it might be more meaningful 
to fix the two arbitrary functions of time ^(t, 0) and ^^{t,0), by requiring that ^(t, rc) and 
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0.2 0.4 0.6 0.8 1 
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Figure 3: The quantity av/{Hor) in Region 1, plotted using parameter values from Table [T] Since 
(Hor^^/c) ~ 0.001, the peculiar velocity av remains small. 



C^itj^c) vanish. This would be in line with the shell r = Tc expanding like the flat EdS 
background. However, this complicates some of the expressions we evaluate, and does not 
change the order of magnitude of any of the final results. Hence we will continue to assume 
that the transformation functions ^ and vanish at r = rather than aX r = Tc- See also the 
end of Section 3.2.] 

The peculiar velocity can be explicitly calculated to be 



a{t)i{t,r) = (Hor) 



1/2 



smu 



1 



cos u 



(42) 



(1 — cosn) 3{u — smu + B 
where the various functions are defined in Eqns. (jlSp . and we have defined the constant B by 



B 



2Hoti 



3/2 



sm Uj 



(43) 



In the rest of this section we will use the parameter values listed in Table [TJ In Fig. [3] we 
have plotted av/{HQr) in Region 1. We see that this dimensionless quantity remains of order 
^ 1 throughout the evolution. For our choice of {H^r^/c) ~ 10~^, which corresponds to an 
overdensity spanning a few Mpc today, the peculiar velocity is of order ~ 10"'^ in Region 1. 
Knowing ^ and we can easily determine and A to be 



7 - (^or)^ 



- -avRH 
2 



(1 



cosn 



smn 



3 \3 (u — sinn + B)"^ [u — sinu + B) 



±p (i/or)Vi(^), 



(44) 



1 f 



^=7 - (^or 



1(2- 



smn 



(1 



cosn 



sin u + B 



sinn + B)"^ 



cos n 



cos n 



{H^rfUu) 



(45) 



where av is given in Eqn. (j42p . and the last equalities in the equations define the functions 
/i(n(t)) and /2(n(t)) respectively. In Fig. |l]we have plotted the functions fi{u{t)) and /2(n(t)). 
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Figure 4: The time dependence of ttie functions ip and A in Region 1, using parameter values from 
Table [TJ Panel (a) shows the evolution from ti to , and panel (b) shows the early time evolution upto 
t = 0.0075^0- The dotted lines correspond to fi{u{t)) and the solid lines to f2{u{t)) (see Eqns. (jH]) and 
P5|) for definitions). We see that at early times, A 2± —^|;. Since {HQr^/c) ~ 0.001, and A remain 
small. 



We see that they remain of order ^ 1 for most of the evolution, and hence for an overdensity 
spanning a few Mpc, and A are of order ~ 10~^ in Region 1. By expanding fi{u) and f2{u) in 
the parameter u around its initial value Ui, one can show that at early times one has A ~ — ^, 
as expected in the Hnear theory. Fig. |4b] shows this behaviour. However, the relative difference 
between A and — grows quickly and (for the parameters given in Table [T]) becomes of the 
order 10~^, by t ~ 0.02to- At the end of Section 3.2 we show that the large relative difference 
between A and —if) (of the order ~ 50% at late times), is largely due to our choice of setting 
these functions to zero at the origin. [We note that a difference between these two functions is 
not in principle surprising since the density contrast in Region 1 grows to order ~ 1 by t ~ 0.3to 
(see Fig. [Tal) . and departures from the linear theory are expected to become significant even 
before this. See, however. Section 4 below. ] 



3.2 The transformation in Regions 2 and 3 

The calculation in Regions 2 and 3 proceeds in a similar fashion as above, but in this case the 
integrals involved cannot be computed analytically. We will therefore display the expressions 
we obtain for i), ^, and ^ and plot the results of numerically computing and A from these 
quantities. 

• Region 2 (r* < r < r^): 

For numerical calculations we found it convenient to define the functions 



Bi{r) 
B2{t,r) 
B3{t,r) 



2 a,- 



IH^ rV 



re 

iTe, 

£3/2 



2 ttj (1 + e) (1 — cos a)2 
Hi sin a 

2^3/2 



Hi[t — ti) smaT:; — ] — '- + 



4^1/2 



sm a 



(l+e) (1 + e) 



sm a,- 



(1 



cos a 



l + e 



(46a) 
(46b) 
(46c) 
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Figure 5: The peculiar velocity av/c in Regions 2 and 3 using parameter values from Table [TJ 



We then have 



v = ^[Bsit,r)-H] , 
a 

C\t,r)=C\t,n) + a{t) [ vit,r)R'{t,f)df, 



(47) 

(48) 
(49) 



where C^{t,r^:) is computed from Eqn. (j4ip at r = r*. ^ and A must now be computed 
using Eqns. (f34]l and (p8l) respectively. We have again used the NIntegrate routine of 
Mathematica. 

Region 3 {vc < r < r^) : 

The analysis is very similar to that in Region 2. We define the functions 



D2{t,r) 
and find 



2 ai 



I Hi rV 



re 

1 + e 



,|3/2 



2 aj (1 + e) (cosh 7/ - 1)2 

i^isinh?? 2|ep/2 
(cosh ry — 1)2 1 + e ' 



N . , (3 + e) 4|e!i/2 g^j^j^^ 
i^i t - smh r,\-—l + -— 
(1 + e) (1 + ey smhj^j 



e(t,r) =e(t,r,)+ r {Di{f) + D2{t,f))dr, 

J Tc 



R 



v = -[D3{t,r)-H] , 
a 



e°(t,r) = e°(t,re) + a(t) / i)(i, f)i2'(i, f)df 



(50a) 

(50b) 
(50c) 

(51) 

(52) 
(53) 



where C(t, r^) and C'^(t, r^) are obtained from Eqns. ()47p and (j49p respectively, evaluated 
in the limit r ^ rr. 
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Figure 6: The metric functions ■(/'(^, J') and A(t^ r) in Regions 2 and 3 using parameter values from Table 
[TJ The time axis begins at t = 50ti. 

In Fig. [5l we have plotted the velocity av/c in Regions 2 and 3 for a range of time. It can be 
shown that at the order of approximation we are working at, av changes sign at r = Tc. [ Recall 
e{rc) = and hence this shell expands exactly like the EdS background. The metric in the 
(t, f) coordinates will not be exactly EdS at r = rc, due to our unusual choice of normalisation 
for ^ and ^'^ at r = 0. This does not pose any problem for our conclusions.] 

In Fig. [6] we plot and A. We see that these functions are well behaved and remain small 
for the entire region of interest (in space and time). Hence the perturbed FLRW picture is 
indeed valid for this system, even though each region by itself appears to be very different from 
FLRW in the synchronous coordinates comoving with the matter. Due to numerical difficulties 
close to the initial time t = ti, we have plotted the time axis starting from t = SOtj. 

Note that the magnitude of ip and A is sensitive to the overall size of the region, determined 
by the value of R{t,ry). For our parameter choices given in Table [H the size of the region at 
the present epoch is ~ 33Mpc, which is a typical size for observed voids. The dependence is 
roughly (HR)^, and hence a void which is about 10 times larger in length scale than the above 
value, would have metric functions about 100 times larger. 

We end this subsection by noting two points. Firstly, we have seen that with our choices 
for certain arbitrary functions of time in Eqns. (j40p and (|4ip. the functions A and build 
up a relative difference of ~ 40-50% by present epoch (see Figs. (Ha|l and ([6])). However, it can 
be shown that if one redefines these functions so that they vanish at r = Vc rather than at the 
origin (in line with the discussion below Eqn. (I4ip ). then the relative difference between the 
functions reduces to ~ 5% at present epoch in Regions 2 and 3. This is demonstrated in Fig. 
[7] for two representative values of the coordinate r. In Region 1 the the relative difference is 
maximum at the center r = 0, being about 30% at present epoch (which is obvious since with 
the redefined functions, the relative difference at the center is simply the relative difference of 
the original functions at r = Vc). The origin of the large relative difference seen in, e.g. Fig. 
l4al is therefore not completely clear, and may even be unphysical. This issue deserves more 
careful consideration, and we hope to return to it in future work (see also Section 4 below). In 
any case, we emphasize that our main result is that the magnitude of the functions A and ip 
themselves, is very small. This brings us to the second issue. 

It is known that simply having a metric of the form (j29p with only the magnitude of the 
perturbations being small, is not enough to guarantee consistency with Einstein's equations 
written as a perturbation series; additional constraints on the derivatives of these functions 
must be satisfied. These constraints, given in e.g. Ref. [8j, take the form (for the metric (I29p 
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Figure 7: The relative difference between the metric functions A and —ip, defined as \A/tp\ — 1, plotted 
as a function of time at r = (l/2)(r, + r^) (panel (a)) and at r = (l/2)(rc + r^) (panel (b)). The solid 
lines show the difference for the definition of the functions used in this paper. The dotted lines show the 
difference when the functions are redefined to vanish on r = Vc instead of vanishing at the origin. 



with ip = —A), 

2 

«C V"iV„i , (v"iV„i) ^ < V^i) V^VfsA , (54) 

where a,f3 = 1,2,3, and Vq, is the spatial covariant derivative associated with the flat 3-space 
metric. On dimensional grounds, treating A ~ {HR)^ <C 1, 5^ ~ if and V ~ aR~^, it is easy to 
see that these constraints will be satisfied by our solution. This should also be expected since 
we started from an exact solution of the Einstein equations and performed a self-consistent 
coordinate transformation. 

3.3 The magnitude of the backreaction 

One can now legitimately ask the question, "How large is the effect of the small metric inho- 
mogeneities?" Naively, one would argue that small inhomogeneities must lead to small effects. 
Indeed, the question of the magnitude of the backreaction in the Newtonianly perturbed FLRW 
setting has been investigated by Behrend, et al. [TO] in the linear and quasilinear regimes, and 
they find that corrections to the FLRW equations remain at the level of one part in 10^. The 
effect of perturbative metric inhomogeneities on observables such as the luminosity distance 
evaluated using the perturbed FLRW metric with a fixed background, has also been studied [11] 
and shown to be small. However, what we are dealing with is a situation in which the matter 
perturbations are completely nonlinear, and it is not a priori clear that the same arguments 
would carry through. Indeed, we saw in section 2 that the deceleration parameter q deviated 
from its EdS value by about ~ 10%. Here we give an argument based on dimensional consider- 
ations supplemented with realistic numbers, which will show that this effect is scale dependent, 
and is not expected to be present if a sufficiently large averaging scale is chosen. 

In the following we will work at the present epoch to- Consider a model situation similar to 
the one we have been considering so far, such that at present epoch the physical extent of the 
overdense region is R^, and that of the underdense is Ry. For order of magnitude estimates, 
we assume that in the perturbed FLRW metric ()29p (which is valid for this system provided 
HoRy <C 1), ^ ~ —V'- Also assume that the density contrast in the overdense region is 
and that in the underdense region is 5yo, where we take (5*o and 5yo to be constant in space, 
which is fine for an order of magnitude estimate. The backreaction in the Buchert approach 
contains, among other terms, the spatial average of the quantity V^A which appears in the 



dA 
'dt 
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spatial curvature [QIIU], where is the Laplacian operator for the flat 3-space metric. The 
spatial curvature has the structure 



7^~ ^ ^ 



(#i)V2i + (#2)iv2i + (#3)(Vi)2 , (55) 



where 7^1, #21 #3 are constants whose values are irrelevant for this order of magnitude argument. 
Due to the Einstein equations in the small scale Newtonian approximation, the leading order 
effect in the nonlinear regime, comes from V'^A which satisfies 

V^i ~ I ^0^*0 , overdense region , 



Hq6.uo 1 underdense region . 

Consider the situation when, at present epoch, i?* ~ 6Mpc, ~ 30Mpc, 5*o ~ 10^ and 
^vo ^ —0.9. These are typical numbers for clusters of galaxies and voids. It is straightforward 
to now show that the spatial average of V'^A over a domain comprising the overdense and 
underdense region, works out to be 

~ -O.IH^ . (57) 

It would appear therefore, that this spatial average of V'^A (which is usually neglected) thus 
turns out to be a significant contributor to the backreaction. (In fact it is the most significant 
contributor, since the other terms are clearly of at least one higher order in the small quantity 
{HqRv)'^, for such a model.) 

As we now argue, however, the above effect can be deceptive, and is really scale depen- 
dent. Let the initial density contrasts in the overdense and underdense regions be and 5yi 
respectively, so that (5*i, \5m\ <^ 1. If M^i, M^, and My are the masses at initial time and 
today, in the overdense and underdense region respectively, and pi and po are the values of the 
background density at initial time and today, then at initial time 

M,i ^ Pi{aiR,f = poRl , M,„^poRl, (58) 

and at present time, 

M^ = Po{l + 6,)Rl> M,i, My = po{l + 5y)Rl < M^i . (59) 

We now make the crucial observation that if the averaging scale is large enough, and we are 
counting several such "pairs" of overdense and underdense regions, then the mass ejected from 
the underdense region must have all gone into the overdense region. It is then easy to show, 
that 

d,Rl « -5yRl , (60) 

which means that, just like in the linear theory, the average of V^A is expected to be negligible 
on such a scale. In the real universe, we do expect that the averaging scale must be at least 
of the order of the homogeneity scale, and on such a scale we will be sampling several pairs 
of overdense and underdense regions. The only cumulative effects that may arise with such 
a choice of scale are from terms such as (VA)^, which as we mentioned earlier, are of one 
higher order in the perturbation and will give effects of the size ~ HQ^H^Ry)^ <C Hq. (For a 
demonstration of the scale dependence of the effect, see e.g. the work of Li and Schwarz |14j.) 

In the light of this discussion, our results of section 2 (figures (pa|) - (pc]) ) can be understood 
as arising due to the fact that we are not averaging on a scale large enough to sample many 
overdense and underdense pairs. In the language of the above argument, all the mass ejected 
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from the underdense region (in coordinates which are comoving with the background), is not 
going into the overdense region which is part of the averaging domain. Some of the matter is 
escaping outside the averaging domain and is left unaccounted for. This explains an imbalance 
as suggested by Eqn. ([57|) . 

4 Discussion 

In this paper we have addressed the question of whether or not the assumptions involved in using 
a perturbed FLRW framework to describe the metric of the present Universe, break down during 
epochs of structure formation, and whether averaging over inhomogeneities in such a context 
can yield effects large enough to lead to an accelerating scale factor. We have demonstrated 
two things. First, it is very essential to impose initial conditions properly matched across any 
boundaries present in the model. Neglecting to do so (for example by ignoring certain regions 
such as Region 2 of our model) can lead to effects such as an accelerating effective scale factor 
amod- Also, as we saw in Fig. I2dl it is in principle possible to have very large deviations from 
FLRW-like conditions after averaging, but this is at the cost of initial conditions which are 
unrealistic when confronted with the CMB data. 

Secondly, we have shown that, provided peculiar velocities of matter are small (which is a 
reasonable assumption on observational grounds), one can explicitly find coordinates in which 
the metric of the Universe is of the perturbed FLRW form (j29p . This by itself is not a surprising 
result, since one can always do this locally (our Region 1 for example is a closed FLRW sub- 
Universe by itself, without any coordinate transformations). What is important is the fact that 
the background FLRW used for comparison is taken to be uniquely and globally defined, so that 
every local region involved in structure formation can be compared to this background in an 
unambiguous fashion. The existence of this global FLRW model is not simply an assumption 
made at the present epoch, it is justified by the fact that density fluctuations in the past were 
very small. 

An issue worth highlighting again is that the time coordinates involved in the transformation 
from the synchronous coordinates to the perturbed FLRW form, are linked by an infinitesimal 
transformation (after controlling some degrees of freedom by setting arbitrary functions of time 
to zero), as expected with weak gravitational fields. While our model is perhaps not the most 
realistic depiction of observed voids, we expect this feature to hold in more realistic models also. 
At first glance this may appear to differ from Wiltshire's claim [3] that the clocks of observers 
in bound regions face a significant calibration when compared with those inside large voids. 
However the issue is rather subtle, since Wiltshire's arguments involve a different choice of time 
slicing than the standard constant t hypersurfaces used here, and his model shows cumulative 
effects over large times. It is possible that these cumulative effects are related to the relative 
difference between the metric functions which we discussed in Section 3, but a straightforward 
comparison is not possible at this stage. In this context, it becomes essential to understand 
the origins of those differences, and moreover the fact that the magnitude of the difference is 
sensitive to boundary conditions. 

A few remarks concerning further tests of the backreaction argument, with more realistic 
models of structure formation which should account for pressure. We have not modelled voids 
surrounded by over densities, since such models are plagued by shell-crossings in the absence 
of pressure. (Our current model also faces this difficulty beyond region 3.) But qualitatively, 
we expect our results to hold even in models where a suitable pressure term takes care of shell 
crossings and leads to stable structures surrounding expanding voids, since we expect peculiar 
velocities to remain small in this situation as well. As long as the sizes of individual voids 
are not an appreciable fraction of the Hubble scale, one expects that here as well, it will be 
possible to find coordinate transformations like the one we have presented. We will return to 
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this problem in future work. 

Finally, we comment on an issue which deserves careful consideration. Recall that in the 
Buchert scheme of averaging perturbative inhomogeneities, which has been employed by several 
authors [3l[ini[l3], there are two scale factors - the background scale factor a{t) which satisfies 
the usual FLRW equations, and the effective scale factor aoit) which satisfies the Buchert 
equations [6] . It is our belief that this situation has an inherent ambiguity - which of the two 
scale factors is relevant for observations? In our opinion, a more consistent way of proceeding 
would be to employ the fully covariant averaging scheme due to Zalaletdinov [7], developed 
further in the cosmological context by [21]. In this scheme, there would be only one scale 
factor, unambiguously associated with a background metric, which would satisfy the corrected 
Einstein equations. The fact that backreaction is negligible in the linear regime should allow 
the problem to be well-posed at say the last scattering epoch. While we have estimated the 
effects of averaging inhomogeneities to be negligible on large enough scales in Buchert 's scheme, 
this was done in a perturbation theory around the fixed Einstein-deSitter background. It is not 
clear that this result will continue to hold even at late times when the background itself is 
"evolving", namely, being corrected by the (small but cumulative) effects of inhomogeneities. 
Note that if these effects build up at late times, then in this scheme with only one scale factor, 
the evolution of the perturbations will also be modified, which may lead to some interesting 
effects. We will investigate these issues in the near future. 

In conclusion, we wish to suggest that it has not yet been conclusively established that 
cosmological backreaction becomes large enough during the structure formation phase to cause 
an acceleration of the scale factor. On the contrary, the calculation presented in the present 
paper suggests that the backreaction remains small during this phase, and that the intuitive 
picture concerning weak gravitational fields is actually realized in the form of a perturbed FLRW 
metric. Self-consistently accounting for the effects that the backreaction would implicitly have 
on the evolution of the perturbations (via the scale factor), may lead to interesting effects which 
remain to be explored. 
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